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ABSTRACT 


A  scheme  for  automatically  integrating  the  heat  equation  is  discussed. 
A  program  based  on  this  scheme  is  explained,  and  the  results  of  integrating  the 
heat  equation  for  a  number  of  different  initial  conditions  are  presented  and 
interpreted. 
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1 .   Introduction 

The  method  of  lines  is  a  technique  for  numerically  solving  partial 
differential  equations  "by  replacing  the  partial  derivatives  of  all  but  one 
independent  variable  by  finite  difference  approximations.   This  results  in 
a  system  of  ordinary  differential  equations  in  the  remaining  independent 
variable  which  can  then  be  solved  by  using  a  standard  numerical  integration 
scheme.   Such  a  technique  is  very  flexible,  since  it  can  be  used  with  almost 
any  type  of  finite  difference  approximations  and  with  almost  any  standard 
method  for  numerically  solving  ordinary  differential  equations.   This  suggests 
that  the  method  of  lines  might  be  well-suited  for  use  in  an  automatic  program 
for  controlling  the  errors  in  the  solution  of  partial  differential  equations. 
By  employing  the  method  of  lines  with  a  fixed  mesh,  an  automatic  program 
could  vary  one  or  more  of  the  following  as  it  deems  it  necessary  in  order  to 
control  the  error  in  accordance  with  any  particular  set  of  criteria: 

1.  The  size  of  the  integration  step  in  the  remaining 
independent  variable. 

2.  The  numerical  integration  scheme  used. 

3.  The  order  of  the  finite  difference  approximations 
used. 

Automatic  programs  for  ordinary  differential  equations  attempt  to 
control  the  global  error  by  bounding  the  local  error  introduced  at  each 
integration  step.   Since  the  theory  behind  this  procedure  is  not  yet  fully 
developed,  it  is  unreasonable  to  expect  the  theory  for  automatically  solving 
partial  differential  equations  to  be  any  further  advanced.   However,  in  the 
following  section,  it  is  shown  that  by  using  the  method  of  lines  and  making 
several  reasonable  assumptions,  it  is  possible  to  employ  an  heuristic  scheme 


to  control  the  error  automat ically  in  the  solution  of  the  heat  equation. 
Section  3  describes  a  program  that  was  written  to  test  experimentally  the 
validity  of  the  analysis  in  section  one  and  describes  the  numerical  tests 
run  with  the  program.   The  final  section  discusses  the  results  of  these 
tests . 


2.   Estimating  the  Global  Error- 


Suppose  that  we  wish  to  solve  numerically  the  heat  equation 

u+  =  Au  (1) 

t     xx 

using  the  method  of  lines.   In  order  to  use  this  technique,  at  t  =  0  we 

choose  a  fixed  set  of  equispaced  points  {x.  }  at  which  we  seek  approximations 

T 
to  the  solution  u(t)  =  (u(t,  x.))   at  the  times  {t.}.   The  t.'s  are  functions 

of  the  local  error  tolerance  e,  since  the  integration  stepsize  must  be 

controlled  in  order  to  keep  the  local  error  below  e.  We  approximate 

m 

u   =  (u  (x.  ,  t))  "by  the  centered  difference  operator  A(T,  e).   Suppose 

that  A(t,  e)  is  controlled  such  that  the  following  conditions  always  hold: 

i.   A(t,  e)  does  not  change  on  [t.,  t  .,-,)•   It  may  change  only 

J    J  •*- 

at  the  points  {t.}.   The  elements  {a.  .(t,  e)}  of  A  are  thus 
piecewise  constant, 
ii.   A(t,  e)  is  symmetric  and  negative  definite  for  all  t  >  0. 
iii.   The  error  p_(t,  e)  =  u   -  Au  in  approximating  the  partial 
derivatives  by  finite  differences  is  controlled  at  t.  such 

J 

that  AAt. I  |a(t . ,  e)  I  I   <  e/2  where  At .  =  t .  n  -  t . .   In 
J   —  J      2  j    j+1    j 

practice,  this  is  too  much  to  expect.   At  the  end  of  this 
section,  we  discuss  the  effect  of  relaxing  this  requirement 
to  make  it  more  realistic. 
Let  A.  denote  A(t.,  e).   In  the  event  that  A( t . ,  e)  4   A(t    ,  e)  define 

J  J  J  J   -*- 

A(t"    e)  =  A(t   e)  =  A 

Let  v  denote  the  exact  solution  to 

v.  =  XA(t,  e)v  (2) 

— t  — 

We  obtain  approximations  w.  to  v.  =  v(t.,  e)  by  numerical  integration  using 

tJ  J  J 

the  trapezoidal  rule 


'  *i+i  =  v-i  +  1/2  XitJ!A(tj^  +  "Vi'^n1 

=  w     +  1/2  XAt  A(t    )[w     +  v        ]  (3) 

J  J  J  J  J       •*- 

Define   t.   to  be  the  local  truncation  error  of  the  trapezoidal  method  in   (3). 
-J 

That   is,   if  w,   =  v.,   then  t.    =  v!   _    -  w...      Assume   for  the  moment  that  we 
-J       -J  -J        -J+l       -J+l 

can   control  the  stepsizes   {At.}   such  that   for  all  j    £  0  ,        x.    L.  =   e/2.      Later, 

J  '  '  -J  '  '  2 

we  will  examine  the  effect  of  slightly  weakening  this  assumption.   Let 

z  =  v  -  ¥  be  the  difference  between  the  exact  and  numerical  solutions  to  (2). 

We  will  now  derive  a  bound  for   z   L. ,  then  address  ourselves  to  the  question 

1  '— n'  '  2 

of  bounding      |u  -  vj  |     ,   and  lastly   combine   the   result  of  the  two   analyses   to 

arrive   at   an   estimate  of  the  global   error    |  |u  -  w|  |      for  the  program. 

In  solving   (2)   by  the  trapezoidal   rule,   each  time  a  step   is   taken, 

the   current  global  error  is  multiplied  by  exp(At.XA.),    and  a  new  local 

J      J 

truncation   error  is    added.      We   thus  have 

z.,_,   =  exp(At .XA. )z.  +   t.  (h) 

-J+l     .  J      J  -J       "J 

Assuming  that  w.  -  V-.  then   z     =   0,   and  applying   (k)  ,  we  have 
"0        — 0  - "0        — 

n-1       n-1 
z     =    i1-      i  •    TLi    exp(At  .XA.  )]t. 

Since  the   A. ' s   are  symmetric  negative   definite,      I exp(At .XA. ) I  I _   $   1.      Using 
l  '  '  J      l       2 

this   fact    and  the   assumption  that      It.IL   =   e/2,  we  obtain 

N^ll2  «  °Joe/2"ne/2  (5) 

We  now  examine      |u  -   vj  |Q.      Using  the   above   definition  of  o_(t,   e), 

we   can  rewrite    (l)    as 

u.    =   XA(t,   e)u  +  Xa(t,   e)  (6) 

— t  —         — 

Subtracting   (2)    and  integrating,  we  have 

t  t 

|u(t    )   -  v(t   )|L   $      |/  n   exp[X/  n  A(s',    e)ds']Xo(s,   e)ds|L 
—     nn^Q  s  —  d 


We   again  use   the   fact   that  A(s',   e)   being  symmetric  negative    definite   implies 
t 

||exp[X/  nA(s',    e)ds']||0   *   x   to  Set 

s  c- 

t  t 

|u(t    )    -  v(t    )|  L    $  /  n    I  |exp[X/  n  A(s»,    e)ds']||l |Xa(s,   e)  |  Lds 

—    n—     n^Q  s  d       —  d 

t 
$   /  n||Xa(s,    e)||_ds 

0  d 

We  approximate  /   J         |  Xo_(  s  ,    e)  |  |    ds  by  At .  X  |  |o_(t .  ,    e)||     .      Assuming  that   the 

t  ^  ^ 

J 

error  in  this    approximation  is    small   enough  to  "be  neglected,  we   can  use  the 

assumption  that    XAt.||o_(t,    e )  |  |       <:   e/2   to    get 

J  <- 

n-1     tA+1 
u(t    )   -   v(t        I      <   e/2    .1.    f  J        ds/At.   =  ne/2  (T) 

'-     n         —     n      '2  J=0   t  J 

J 

Finally,   we   combine   (5)    and  (7)    to   get   a  bound  for      |u  -  w|  |     .      Using  the 


triangle   inequality,   we  get 


|u  -  w|  |      $  ne  (8) 


Since  we  are  using  the  trapezoidal  rule  to   integrate   (2)   numerically,    it   is 

reasonable  to  expect  that      Ix-llp  =   0(AtT>).      But,    since      |  t_.  |  |      =   e/2,    it 

3  1/3 

follows   that  At.   =  0(e)    and  hence,   that  At.   =   0(e        ).      If  we  now  assume  that 

J  J 

-1/3 
At.  =  0(l/n)  (as  it  ought  to  be),  then  we  have  n  =  0(e     ).   Substituting 

J 

this  relation  in  (8),  we  get 


|u  -  w||2  $  ke2/3  (9) 


for  some   constant  k. 


In  practice,  we   cannot   expect   a  program  to   actually   control  o_(t.,    e) 

J 

and  x_.  as  strictly  as  was  assumed  in  the  above  analysis.   We  will  now  look 

at  the  consequences  of  reducing  these  expectations  to  a  realizable  level. 

Since  u  is  smooth  and  A(t,  e)  is  piecewise  smooth  in  t,  we  can  express  the 

truncation  error  of  the  trapezoidal  method  at  t .  ,  as  a  power  series  in  At. 

J+l  J 


00 

x.  =    e    x.  .At: 

-j  1=3   -tj       0 

Instead  of  requiring  that    ||t_JIo  =   £/25  we   require   that    I|t     At.||      $   e/2. 

In  practice,  we   do  not  know   t_       exactly,   but  we  will   assume  that  we   can 

estimate  it   accurately  enough  so  that  the   following  conclusions    are   not 

altered.      Since      |  t_    A"t .  |  L  $   e/2,  we  may  infer  that  At.   =  0(e   1/3),    and 

-30       J  u 

hence 

|  1  t_  |  |      =    |  |t     At^|  I      +   o(At^)    *   e/2  +   o(e)  (10) 

The    (at  least)   piecewise   smoothness  of  u  and  A(t,   e)    allows  us  to 

00  "1 

similarly  write  a(t,  e)  on  [t.,  t.,n)  as  a(t,  e)  =  .1   a.(t,  e)Ax  .   We 

-  j    o+l     -        T=p.  -l 

will  no  longer  purport  to  control  o_   such  that  XAt .  |  |£.(t.,  e)||   <  e/2. 

J  J  *- 

Instead,  we  only  require  that 

Pi,  , 
XAt.   o   (t.,  e)Ax  J  L  $  e/2  (ll) 

0  '-p.   0  2 

J 

Once  again,  in  practice,  we  only  have  an  estimate  of  a   (t.,  e),  but  we  will 

assume  that  the  estimate  is  sufficiently  accurate  for  the  discussion  to 

1/3  2/3pl 

remain  unaltered.   Since  At.  =  0(e    ) ,  it  follows  from  (ll)  that  Ax  =  0(e    J ) 

0 

and  hence, 

1+2/3P, 

At.  |  |a(t,  e)|  |  S  e/2  +  o(e      J  )  (12) 

J 

We  can  therefore,  still  expect  the  results  of  the  previous  analysis  to  hold, 
even  when  the  requirements  are  less  stringent. 


3.   Programming  Application 

A  version  of  the  numerical  procedure  analyzed  in  the  previous  section 
was  implemented  in  a  FORTRAN  IV  program.   This  section  "begins  by  giving  some 
of  the  particulars  of  its  operation.   For  a  more  detailed  account  of  its 
functioning,  see  the  block  diagram  which  has  been  included  in  the  appendix. 
The  program  was  written  specifically  to  solve  the  heat  equation  for  various 
initial  conditions.   It  approximates  the  second  spatial  derivatives  using 
centered  differences  and  uses  the  trapezoidal  rule  to  integrate  with  respect 
to  time.   It  was  numerically  verified  that  all  of  the  centered  difference 

formulas  used  give  rise  to  symmetric  matrices  A(  t  ,  e)  which  are  negative 

J 

definite. 

Prior  to  attempting  an  integration  step,  an  estimate  of  the  local 
spatial  truncation  error  is  made  using  centered  differences  on  the  current 
numerical  solution  to  estimate  the  magnitude  of  the  higher  order  derivative 
which  appears  in  the  leading  term  of  the  Taylor  expansion  for  the  error. 
If  the  estimated  error  exceeds  e/2 ,  the  current  centered  difference  formula 
for  the  second  derivatives  is  replaced  by  another  centered  difference 
formula  which  is  two  orders  of  accuracy  in  Ax  higher.   The  local  spatial 
truncation  error  for  the  new  formula  is  estimated,  and  if  its  magnitude 
still  exceeds  e/2,  the  process  repeats  with  the  centered  difference  formula 
that  is  higher  in  order  by  yet  two  more  powers  of  Ax.   This  continues  until 
either  a  formula  is  found  for  which  the  estimated  error  is  less  than  e/2  in 
magnitude  or  the  order  of  the  formula  reaches  a  limit  specified  in  the 
initialization  section  of  the  program. 

If,  on  the  other  hand,  the  estimated  spatial  truncation  error  is  much 
less  than  e/2,  an  estimate  is  made  of  what  the  spatial  truncation  error  would 


8 

be   if  the    current   centered  difference   formula  were   replaced  by  a  centered 
difference   formula  whose   accuracy  is   two   orders   lower.      If  the   estimated 
error   for  this   lower  order  formula  still  is   less  than   e/2  in  magnitude,   the 
lower  order   formula  replaces   the   current   formula. 

An   integration   step  in  t   is  then  taken,    and  the  temporal  truncation 
error  is   estimated  using   step-doubling.      This   procedure   involves  taking  a 
step  of  size  At   and  comparing  the  result  with  the   solution  obtained  by 
taking  two  steps   of  size  At/2.      If  the   estimated  error  exceeds   e/2,   At  is 
adjusted  downward,    and  the  integration  step  is   reattempted  with  the  adjusted 
At.      If,   instead  the   estimated  error  is   less   than   e/2,   the   step  is    accepted, 
and  At  may  possibly  be   increased  for  the  next   step. 

The   program  was   used  to  solve  the  heat  equation 

u,    =   Xu 
t  xx 

with  the  diffusivity  constant  X  =  1/10  subject  to  the  initial  condition 

u(x,  0)  =  f(x) 
We  assumed  f(x)  was  periodic  with  period  2  and  antisymmetric  with  respect  to 
x  =  0,  ±1,  ±2,  ...  .  Hence,  we  can  restrict  our  attention  to  the  interval 
[0,  l].   This  reduces  the  run  time  and  also  permits  us  to  use  centered 
differences  at  all  mesh  points  in  [0,l],   Under  these  assumptions,  with 
boundary  conditions 

u(0,  t)  =  u(l,  t)  =  0, 

the  solution  is 

oo  2   2 

u(x,   t)    =     Z-   b     exp(-n   tt   Xt)    sin  nnx 

n=l     n 

where  b     =  2  /      f(x)   sin  nirx   dx 

n  0 

The   following  initial   conditions   were   used: 


Case   1.      f(x)   =  sin  ttx, 

b     =   1,        n  =   1, 
n 

0,        n   >   1. 
Case  2.      f(x)   =1        0    <  x   <  1. 
b     =   U/(m0  ,        n   odd, 
0,        n   even. 


Here   u(x,   0)    is    a  square-wave. 

Case   3.      f(x)   -  kx,        0    $  x   $   .25, 

1,        .25    S  x  *   .75, 

h  -  kx,       .75  $  x  $  l. 

b     =      8   (sin  mr/U   +  sin   3nir  A ) ,        n  odd, 
n       — 
nu 

0,        n  even. 

Here   u(x,    0)    is    an   isoceles   trapezoid  on    [0,   l]. 
Case   U.      f(x)   =  x/a,        0   <  x   £   a, 

=   1  -   (x  -   a)/(l  -   a),        a  {  x  $' 1. 

2   2 
b      =   2    sin    (mrp)/[n    tt   p(l   -a)]. 

.  Case  Ua.  p  =   1/2 

Case  Ub.  p  =   7/10 

Case  Uc.  p  =   9/10 

Case  Ud.  p  =  999/1000 

Here  u(x,   0)   on   [0,   l]    is    a  peak  with  its   summit   at   x  =    a.      As    a 
approaches   1,   the  slope  to  the  left   of  p  slowly  decreases,  while 


10 

the  slope  to  the  right  of  y  becomes  steeper  at  an  ever  increasing 
rate.   For  y  close  to  1,  u(x,  0)  is  very  nearly  a  sawtooth  wave. 

_2 
The  program  was  run  for  each  of  the  above  cases  with  e's  of  10   , 

-3        -8  -9 

10   ,  . .  . ,  10   except  for  CASE  1  where  it  was  also  run  for  e  =  10   and 

e  =  10    .   For  all  cases  except  for  CASE  1,  the  integration  was  begun  at 

t  =  1/10  rather  than  at  t  =  0  in  order  to  eliminate  the  difficulties  arising 

from  the  discontinuities  at  zero.   The  initial  number  of  points  in  the 

centered  difference  approximation  to  the  second  derivative  was  three.   The 

initial  settings  for  At  and  the  number  of  internal  mesh  points  were  .01  and 

19  respectively.   Integration  was  halted  as  soon  as  t  exceeded  2. 


11 

k.      Results 

The   results   of  the  tests    described  above   are  summarized  in   the   tables 
which  appear  in   this    section.      An  explanation  of  each   column  that   appears   in 
the   table   is   given  below: 
e_  -  On  each  integration  step,   the  magnitude  of  the   estimated  spatial   and 

temporal  truncation  errors    are  both   required  to  be   less   than  or  equal 

to  e/2. 
#_  OATT.fl  to  RKSTEP  -  Each  time  that   the   subroutine  RKSTEP  is    called,    an 

integration  is  performed  using  the  trapezoidal   rule.      Since  the 

program  uses   step   doubling,    each  time   a  step   is   attempted,   the  number 

of  calls   to   RKSTEP  is    increased  by  three. 
#_  CALLS   times   #_  POINTS  -  Each   time  RKSTEP  is    called,   the   current  number  of 

internal  mesh  points  is   added  to   this   number.      It   thus   provides   a 

rough  measure  of  the  work   done  by  the   program. 
LOG  10  DIFFERENCE   in  WORK  -  The   difference  between  the  log   (base  10)    of  the 

entry   from  the   current   row  of  the   tf  CALLS   TIMES   #  POINTS   column   and 

the  log   (base  10)   of  the   entry   from  the  preceding  row  in  that   column. 

This  may  give   a  better  idea  of  the  rate  at  which   the  work   done  by  this 

program  increases    for  each   factor  of  ten   decrease   in   e. 

INTERNAL  MESH  POINTS  -  The  number  of  equispaced  points    {x. }    (exclusive   of 

x^  =   0   and  x    ,_    =  l)    at  which   the   solution   is   computed.      If  the 
0  n+1 

program  has    difficulty  in  taking  the   first  step  successfully,   it 
may  attempt  to  raise   the  number  of  points  . 
INITIAL  COUPLING  -  The  number  of  points  used  in  the   centered  difference 
approximation  to  the   second  derivative  when  the   first   successful 
step  was   taken. 
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FINAL  COUPLING  -  The  number  of  points  used  in  the  centered  difference 
approximation  to  the  second  derivative  when  the  upper  limit  of 
integration  was  reached. 
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FINAL  At  -  The  stepsize  when  the  upper  limit  of  integration  was  reached. 

MAX  ERROR  EVER  -  The  magnitude  of  the  worst  error  in  any  individual  component 

of  the  numerical  solution  vector  that  was  ever  found.   The  numerical 

solution  was  compared  to  the  exact  solution  on  each  step  that  t 

equalled  or  exceeded  an  integral  multiple  of  2/10. 
MAX  FINAL  ERROR  -  The  magnitude  of  the  worst  error  among  the  components  of  the 

final  numerical  solution  vector. 
LOG  10  DIFFERENCE  IN  MAX  FINAL  ERROR  -  This  column  is  to  MAX  FINAL  ERROR  as 

LOG  10  DIFFERENCE  IN  WORK  is  to  §   CALLS  TIMES  #  POINTS. 

Results  indicate  that  the  idea  of  controlling  the  global  error  through 
separate  control  of  the  local  spatial  and  temporal  discretization  errors  is 
indeed  quite  feasible.   Of  particular  interest  is  the  LOG  10  DIFFERENCE  IN  MAX 

FINAL  ERROR  columns.   These  show  that  for  each  factor  of  ten  decrease  in  e, 

2/3 
the  global  error  decreased  by  a  factor  which  appeared  to  approach  (10) 

This  is  as  predicted  in  the  previous  section. 

The  LOG  10  DIFFERENCE  IN  WORK  columns  show  that  as  e  was  successively 

decreased  by  factors  of  ten,  the  amount  of  work  done  increased  by  factors 

1/3        1/2 
which  tended  to  stabilize  in  the  range  of  (10)    to  (10)    with  the  exact 

number  depending  on  the  initial  condition.   For  CASE  1,  where  it  was  never 

found  necessary  to  increase  the  number  of  internal  mesh  points,  the  number  was 

1/3 
precisely  (10)    .   This  is  what  would  be  expected  with  the  trapezoidal  rule. 

Since  the  local  truncation  error  is  proportional  to  the  cube  of  the  stepsize, 

and  the  program  attempts  to  keep  the  local  errors  on  the  order  of  e,  if  the 

-1/3 
mesh  is  fixed,  the  number  of  integration  steps  performed  should  be  0(e    ). 

Combining  the  relationships  between  e  and  the  global  error  and  between 

e  and  the  amount  of  work  done  by  the  program,  we  can  get  a  relationship 
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between  the  global  error  and  the  work.   Denoting  the  global  error  by  e  and 

the  work  by  q,  we  can  write 

2/3 
limit  e(e)  =  c  e 

E  -*■    0 

limit   q(e)   =    c   e~m  1/3   ^   m  $   1/2 

e  -*■  0 

and  hence, 

limit    e(q)    =    c   qP  -2   £  p   $    -k/3 

e   ->  0  J 

where   c.  ,    c_  ,    and  c      are    constants. 

While   this   may  not   seem  like   a  particularly  good  return   of  accuracy 

for  the   computational  effort   expended,    it    compares    quite   favorably  with  most 

of  the  methods    commonly  used  to  solve  partial   differential   equations.      The 

well-known  Crank-Nicolson  method  which  for  the  heat   equation  gives   rise  to 

the   system 

n+1  n        .  ,_.    ,    /  .2     n         „2     n+lN 

u.        =  u.    +  1/2   Ar( 6     u.    +    6     u.       ) 
i  i  x     l  x     i 

where  u.    =  u(x. ,   t),   r  =   At/Ax    ,   and   6      is  the   standard  3-point   centered 

11  x 

3  2 

difference,  has   a  local  truncation  error  of  0(At     +  At   Ax   ) .      The  Douglas 

formula   is  the  most   accurate   formula  involving  the  same   six  points   to   approx- 
imate u     (x.,   t    ._).      It   generates  the   system 
xx     i        n+1 

n+1   _     n  ,  ,  x2     n    ,     x2     n+lx  1   ,,.2     n        x2     n+lx 

u.         =  u.    +   1/2   Xr[&     u.    +   6     u.       )    +  r-r  X(  6     u.    -   6     u.       ) 
i  i  xi  x     l  12  xi  xi 

for  the   heat   equation  which  results    in   a  local   error-  of  0(At     +  At  Ax   ) .      If, 

as  At  approaches   zero,    r  is   kept   fixed,   the  local   error  is   0(At    ) .      The  global 

2 
error  should  therefore  be  0(At    ).      This  means   that 

limit   e(At)   =   0(At2) 

At  ->  0 

just  as  for  the  trapezoidal  rule. 

In  advancing  the  solution  one  time  step,  Douglas'  method  must  solve  a 
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tridiagonal  system.   The  amount  of  work  necessary  to  solve  such  a  system  is 
proportional  to  the  number  of  equations  in  the  system  and  hence  to  the  number 
of  internal  mesh  points.  But  since  the  number  of  internal  mesh  points  is 
always  within  one  of  Ax   ,  we  may  take  the  work  proportional  to  Ax   .   The 
number  of  steps  that  must  be  taken  and  hence,  the  number  of  times  the  system 
must  be  solved  is  proportional  to  At   .   Hence,  for  Douglas'  method,  Ax   At 
may  be  taken  as  a  measure  of  work  analogous  to  the  measure  #  CALLS  TIMES  # 
POINTS  described  above  for  the  program.   If,  as  before,  r  is  assumed  to  remain 
fixed  as  At  varies,  for  the  Douglas  method  we  have 

limit  q(At)  =  0(Ax_1  At"1)  =  0(At~3'2) 

At  ■*■   0 


and  hence, 


limit  e(q)  =  0(q  /3) 


At  +   0 
which  is  no  better  than  the  worst  case  observed  for  the  program. 
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Print  computed 
solution  in 
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